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lO ! Abstract 

Q . An algorithm to estimate motion from satellite imagery is presented. 

O . Dense displacement fields are computed from time-separated images of of 



significant convective activity using a Bayesian formulation of the motion 
estimation problem. Ordinarily this motion estimation problem is ill-posed; 
there are far too many degrees of freedom than necessary to represent the 



/\ • motion. Therefore, some form of regularization becomes necessary and by 



imposing smoothness and non-divergence as desirable properties of the es- 
timated displacement vector field, excellent solutions are obtained. Our ap- 
proach provides a marked improvement over other methods in conventional 
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1 INTRODUCTION 2 

use. In contrast to correlation based approaches, the displacement fields 
produced by our method are dense, spatial consistency of the displacement 
vector field is implicit, and higher-order and small-scale deformations can 
be easily handled. In contrast with optic-flow algorithms, we can produce 
solutions at large separations of mesoscale features between large time-steps 
or where the deformation is rapidly evolving. 

1 Introduction 

Environmental data assimilation is the methodology for combining imperfect model 
predictions with uncertain data in a way that acknowledges their respective un- 
certainties. The proper framework for state estimation includes sequential lfT51l . 
ensemble-based [14J and variational li20l l5l methods. 

The difficulties created by improperly represented error are particularly appar- 
ent in mesoscale meteorological phenomena such as thunderstorms, squall-lines, 
hurricanes, precipitation, and fronts. We are particularly interested in rainfall data- 
assimilation, where rainfall measurements from satellite data, radar data, or in-situ 
measurements are used to condition a rainfall model. Such conditional simula- 
tions are valuable both for producing estimates at the current time (nowcasting), 
as well as for short-term forecasting. 

There are a countless number of models developed to simulate the rainfall 
process. In general, there are two types of models that can deal with spatial and 



1 INTRODUCTION 3 

temporal characteristics of rainfall. The first category is the meteorological model 
or the quantitative precipitation forecasting model. It involves a large, complex 
set of differential equations seeking to represent complete physical processes con- 
trolling rainfall and other weather related variables. Examples of these models in- 
clude the fifth-generation Mesoscale Model (MM5) li3l|4l|T6l|, the step-mountain 
Eta coordinate model dlZlIISj and the Regional Atmospheric Modeling System 
(RAMS) dim, etc. The second type is the spatio temporal stochastic rainfall 
model. It aims to summarize the spatial and temporal characteristics of rainfall by 
a small set of parameters |l6l [TSl [Til HI l22l IZHl- This type of model usually sim- 
ulates the birth and decay of rain-cells and evolve them through space and time 
using simple physical descriptions. Despite significant differences among these 
rainfall models, the concept of propagating rainfall through space and time are 
relatively similar. 

The major ingredient required to advect rainfall is a velocity field. Large 
spatial-scale (synoptic) winds are inappropriate for this purpose for a variety of 
reasons. Ironically, synoptic observations can be sparse to be used directly and 
although synoptic- scale wind analyses produced from them (and models) do pro- 
duce dense spatial estimates, such estimates often do not contain variability at the 
meso-scales of interest. The motion of mesoscale convective activity is a natural 
source for velocimetry. Indeed, there exist products that deduce "winds" by esti- 
mating the motion of temperature, vapor and other fields evolving in time llQlfTOl. 
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In this paper, we present an algorithm for velocimetry from observed motion 
from satellite observations such as GOES, AMSU, TRMM, or radar data such as 
NOWRAD. This algorithm follows from a Bayesian formulation of the motion es- 
timation problem, where a dense displacement field is estimated from two images 
of cloud-top temperature of rain-cells separated in time. Ordinarily, the motion 
estimation problem is ill-posed, because the displacement field has far too many 
degrees of freedom than the motion. Therefore, some form of regularization be- 
comes necessary and by imposing smoothness and non-divergence as desirable 
properties of the estimated displacement vector field solutions can be obtained. 

This approach provides marked improvement over other methods in conven- 
tional use. In contrast to correlation based approaches used for deriving velocity 
from GOES imagery, the displacement fields are dense, quality control is implicit, 
and higher-order and small-scale deformations can be easily handled. In contrast 
with optic-flow algorithms lEDEI, we can produce solutions at large separa- 
tions of mesoscale features between large time-steps or where the deformation is 
rapidly evolving. 

After formulating the motion estimation problem and providing a solution, 
we extend the algorithm using a multi-resolution procedure. The primary ad- 
vantage of a multi-resolution approach is to produce displacement fields quickly. 
The secondary advantage is to structure the estimation homotopically; coarse or 
low-frequency information is used first to produce velocity estimates over which 
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deformation adjustments from finer-scale structures is superposed. The result is a 
powerful algorithm for velocimetry by alignment. As such, it is useful in a variety 
of situations including, for example, (a) estimating winds, (b) estimating transport 
of tracers, (c) Particle Image Velocimetry, (d) Advecting Rainfall models etc. 

2 Related Work 

There are two dominant approaches to computing flow from observations directly. 
The first is correlation-based and the second is based on optic flow. 

In correlation based approaches ifT^ . a region of interest (or patch) is identified 
in the first image and correlated within a search window in the second image. The 
location of the best match is then used to compute a displacement vector. When 
the input image or field is tiled, possibly overlapping, and regions of interest are 
extracted from each tile location, the result is velocimetry at regular intervals and 
is most commonly used for Particle Image Velocimetry (PIV). In certain instances 
it is useful to define interest-points or salient features around which to extract re- 
gions of interest. In particular, if the field has many areas with negligible spatial 
variability, then matches are undefined. As a quality control measure then, match- 
ing is restricted only to those regions of interest that have interesting variability, 
or interest points. 

There are several disadvantages to correlation-based approaches. First, by 
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construction it is assumed that the entire ROI purely translates from one image to 
the other. This is not always the case, but is a reasonable approximation when the 
right length scale can be found. However, when higher-order deformations (shears 
for example) are present, correlation based approaches cannot be expected to work 
well. Second, correlation based approaches assume that a unique match can be 
found in a way that is substantially better than correlation elsewhere. This is only 
true if the features are well-defined and identified. Third, there is no implicit con- 
sistency across regions of interest in correlation-based flow. Neighboring regions 
of interest can and often do match at wildly different and inconsistent locations. 
This calls for a significant overhead in terms of quality control. Fourth, it is not 
clear how the search window size (that is the area over which a region of interest 
is matched in the subsequent frame) is determined. This window size varies both 
in space (as the velocity varies spatially) and time (as velocity varies with time). 
A larger search window portends a larger probability to miss the real target, and 
a smaller search window can lead to false negatives or false positives. Finally, 
where interest points are used as a preprocessing step to correlation, the velocity 
field produced is necessarily sparse, and therefore, leaves hanging the question of 
how to produce dense flow fields. Our proposed algorithm handles all these issues 
in a simple and direct way. 

More closely related to the proposed approach is optic flow lETl \T7\ . This 
method arises from what is known as the brightness constraint equation, which 



2 RELATED WORK 7 

is a statement of conservation of brightness (intensity) mass, expressed by the 
continuity equation evaluated at each pixel or grid node of X. 



dX 

— + q-VX = (1) 

ot 



Here X is the brightness or intensity scalar field and q a displacement vector- 
field. Solutions to the optic flow equation can be formulated using the well-known 
method by ll2Tl . which can be stated as a solution to the following system of 
equations: 

f)X 

(VX)(VXfq=-(VX)— (2) 

The right-hand side is completely determined from a pair of images and the 
coefficient or stiffness matrix on the left-hand side is the second-derivative of the 
auto correlation matrix, also known as the windowed second-moment matrix, or 
Harris interest operator, which is sensitive to "comers" in an image. This formula- 
tion arises directly from a quadratic formulation, which can in turn be synthesized 
from a Bayesian formulation under a Gaussian assumption. Thus, we can write 
that we seek to minimize 

J(q) = ||X(r-q)-y|| (3) 
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Then solve this problem via the Euler-Lagrange equation: 

^^ = VX|r_q(X(r-q)-r) = (4) 

The solution is obtained by linearizing @, that is, 

VX|r-q(X(r)- VX-q-y) = 

VXiVXfq, = -VX{Y-{X{r)) (5) 

There are several disadvantages to this algorithm. First, much like correlation 
with feature detection, equation|5]is evaluated at pixels where the second-moment 
matrix is full-rank, which corresponds to locations where features are present. 
There is no clear way of propagating information obtained at sparse locations to 
locations where direct computation of displacement is not possible due to poor 
conditioning of the second-moment matrix. For the same reason, it cannot han- 
dle tangential flows. The brightness constraint equation can only represent flows 
along brightness streamlines. When tangential motion is present, detected mo- 
tion at extreme ends a moving curve cannot be propagated easily into the interior. 
Our method provides some degree of spatial smoothness common in geophysical 
fluid transport, and uses regularization constraints to propagate flow information 
to nodes where feature strengths are weak. 

Second, the linearization implicit in ^ precludes large displacements; struc- 
tures must be closely overlapping in successive images, which can also be seen 
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from the continuity equation ^. Therefore, this method is very useful for densely 
sampled motion, such as ego-motion resulting from a moving, jittering camera, 
but is not as useful for sparsely sampled flow arising from structures moving in a 
scene. In the latter case, to ameliorate the effects of large expected displacement, 
multi-resolution approaches have been proposed. Even so, much like determining 
the size of the search window in correlation, determining the number of resolu- 
tions is an ad-hoc procedure. Our method can handle large displacements and 
we also propose a multi-resolution approach, but the primary motivation there is 
improved computational speed. 

3 Velocimetry by Field Alignment 

The main approach consists of solving a nonlinear quadratic estimation problem 
for a field of displacements. Solutions to this problem are obtained by regularizing 
the an ill-posed inverse problem. The material presented in this section is derived 
directly from work by Ravela lE^ . and Ravela et al. lE^ . Here we reformulate 
their original formulation to allow only position adjustments. 

To make this framework more explicit it is useful to introduce some nota- 
tion. Let X = X(r) = {X[r[] . . .X[rJ^]} be the first image, written as a vec- 
tor, defined over a spatially discretized computational grid il, and r^ = {r, = 
(xj, Vi)^ , i E Vl}he the position indices. Let q be a vector of displacements, that 
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is q""" = {q. = (Axi, Ayi)'^ ,i E Vl}. Then the notation X(r — q) represents 
displacement of X by q. The displacement field q is real- valued, so X(r — q) 
must be evaluated by interpolation if necessary. It is important to understand that 
this displacement field represents a warping of the underlying grid, whose effect 
is to move structures in the image around, see Figure [TJ 
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Figure 1 : A graphical illustration of field alignment. State vector on a discretized 
grid is moved by deforming its grid (r) by a displacement (q). 



In a probabilistic sense, we may suppose that finding q that has the maximum 
a posteriori probability in the distribution P(q| Af, 3^) is appropriate. Without loss 
of generality, A" is a random variable corresponding to the image or field at a given 
time and y is random variable for a field at a future time. Using Bayes rule we 
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obtain P{Q = c[\X = X,y = Y) oc P{y = Y,X = X|q)P(q). If we make a 
Gaussian assumption of the component densities, we can write: 

p(X, Flq) = i _g-i(y-X(r-q))-i?-l(y-X(r-q)) (g) 

(27r)t |i?|2 

This equation says that the observations separated in time can be related using 
a Gaussian model to the displaced state X(r- q), where X(r) is defined on the 
original grid, and q is a displacement field. We use the linear observation model 
here, and therefore, Y = HX{y ^ q) + r],rj r^ N{0, R).. We should emphasize 
here that the observation vector is fixed. It's elements are always defined from 
the original grid. In fully observed fields, H is an identity matrix, and for many 
applications R, reflecting the noise in the field, can also be modeled as an identity 
matrix. 

^(q) = ^e-"^^^^ (7) 

This equation specifies a displacement prior. This prior is constructed from 
an energy function L(q) which expresses constraints on the displacement field. 
The proposed method for constructing L is drawn from the nature of the expected 
displacement field. Displacements can be represented as smooth flow fields in 
many fluid flows and smoothness naturally leads to a Tikhonov type formulation 
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and, in particular, L(q) is designed as a gradient and a divergence penalty 
term. These constraints, expressed in quadratic form are: 

^(q) = ^ E tr{[Vg;[Vg/} + f E[V ■ g/ (8) 

In EquationlH q^ refers to the f^ grid index and tr is the trace. Equation[8]is 

a weak constraint, weighted by the corresponding weights wi and W2. Note that 

the constant C can be defined to make Equation IT] a proper probability density. In 

particular, define Z(q) = e^^^*^^ and define C = J Z(q)dq. This integral exists 

q 

and converges. 

With these definitions of probabilities, we are in a position to construct an 
objective by evaluating the log probability. We propose a solution using Euler- 
Lagrange equations. Defining p = r — qThese can be written as: 



^ = vx|pi7-i?-(//x(p)-y) + ^ = o (9) 



Using the regularization constraints (|9l) at a node i now becomes: 

WiV'^q. + W2V{y ■ q) + \VXf^\pH^R-^ (h Ix^ (p)! - y)j . = (10) 



Equation [TOl is the field alignment formulation. It introduces a forcing based 
on the residual between the model- and observation-fields. The constraints on the 
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displacement field allow the forcing to propagate to a consistent solution. Equa- 
tion[TOlis also non-linear, and is solved iteratively, as a Poisson equation. During 
each iteration q is computed by holding the forcing term constant. The estimate of 
displacement at each iteration is then used to deform a copy of the original fore- 
cast model-field using bi-cubic interpolation for the next iteration. The process 
is repeated until a small displacement residual is obtained, the misfit with obser- 
vations does not improve, or an iteration limit is reached. Upon convergence, we 

N 

have an aligned image X{p), and a displacement field q = Z) q^^\ for individual 

d=i 

displacements q^'^^ at iterations d = 1 . . . D 

3.1 Multi-resolution Alignment and Velocimetry 

The convergence of solution to the alignment equation is super-linearly dependent 
on the expected displacement between the two fields. Therefore, it is desirable to 
solve it in a coarse-to-fine manner, which serves two principal advantages. The 
first, as the following construction will show, is to substantially speed-up the time 
to alignment because decimated (or coarse-resolution) representations of a pair of 
fields has smaller expected displacement than a pair at finer resolution. 

Second, decimation or resolution reduction also implies that finer structure or 
higher spatial frequencies will be attenuated. This smoothness in the coarsened- 
field intensities directly translates to smoothness in flow-fields using (|9l). Thus, a 
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coarse-to-fine method for alignment can incrementally add velocity contributions 
from higher-frequencies, that is it incrementally incorporates higher-order vari- 
ability in the displacement field. Many of the advantages of a multi-resolution 
approach have been previously explored in the context of visual motion estima- 
tion, including the famous pyramid algorithm and architecture for matching and 
flow and our implementation borrows from this central idea. 
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Figure 2: The multi-resolution algorithm is shown for two-levels and requires five 
steps, labeled (1) through (5). See text for explanation. 



The multi-resolution algorithm is depicted in Figure |2| for two levels. The 
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input images X and Y are decimated to generate coarse resolution images Xi 
and Yi respectively (step 1). Let us suppose that this scaling is by a factor of 
< s < 1 (most commonly s = 0.5). Displacement is computed for this level 
first, and let us call this qi (step 2). This displacement field is downscaled by 
a factor of s, using simple (bicubic) interpolation, to produce a prior estimate of 
displacement at level 0, written qio = s^^qo(s^^r) (step 3). The source image at 
level 0, that is Xq = X is displaced by qio (step 4) and thus X{r — qio) is aligned 
with Yq to produce a displacement estimate qo (step 5). The total displacement 
relating source image X with target field Y is simply qo + qio- Multiple levels 
of resolution can be implemented from this framework recursively. 

4 Example 
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Figure 3: CIMSS Winds derived from GOES data at 2006-04-06-06Z (left) and 
pressure (right). The velocity vectors are sparse and contain significant diver- 
gence. 
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Figure 4: CIMSS Winds derived from GOES data at 2006-04-06-09Z (left) and 

pressure (right). The velocity vectors are sparse and contain significant diver- 
gence. 



The performance of this algorithm is illustrated in a velocimetry computation. 
To compare, we use CIMSS wind-data satellite data ifTDI . depicted in Figure|3J and 
Figure m obtained from CIMSS analysis on 2006-06-04 at 06Z and 09Z respec- 
tively. CIMSS wind-data is shown over the US great plains, and were obtained 
from the 'sounder.' The red dots indicate the original location of the data. The left 
subplot shows wind speed (in degree/hr). The right ones show pressure, and the 
location of raw measurements in red. 

It can be seen in the maps shown in Figure|3and FigurelUthat current method 
to produce winds generate sparse vectors and, further, has substantial divergence. 
Whilst this can be thought of as accurately representing turbulence, in reality these 
vectors are more likely the result of weak quality control. The primary methodol- 
ogy used here is to identify features in an image, extract regions of interest around 
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them and search for them in subsequent frames. This, by definition produces 
sparse velocity estimates (features are sparse), leaving unanswered how to sys- 
tematically incorporate appropriate spatial interpolation functions for the velocity. 
Since regions of interest are essentially treated as being statistically independent, 
mismatches can produce widely varying displacement vectors. Such mis-matches 
can easily occur in correlation based approaches when the features are not distin- 
guishing or substantial deformations occur from one time to another in a region 
of interest. A more detailed discussion is presented in Sectional 

In contrast, our method produces dense flow fields, and quality control is im- 
plicit from regularization constraints. Figure |5ta,b) shows a pair of NOWRAD 
images at 2006-06-0 1-0800Z and 2006-06-0 1-0900Z respectively, and the com- 
puted flow field in Figure |5tc). Similarly, Figure |5td,e,f) show the GOES images 
and velocity from the same time frame over the deep convective rainfall region in 
the Great Plains example. The velocities are in good agreement with CIMSS de- 
rived winds where magnitudes are concerned, but the flow-fields are smooth and 
visual confirmation of the alignment provides convincing evidence that they are 
correct. 
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5 Conclusions 

Our method is a Bayesian perspective of the velocimetry problem. It has several 
distinct advantages: (a) It is useful for a wide range of observation modalities, 
(b) Our approach does not require features to be identified for computing velocity. 
This is a significant advantage because features cannot often be clearly delineated, 
and are by definition sparse, (c) Our approach implicitly uses quality control in 
terms of smoothness, and produces dense flow-fields, (d) our approach can be 
integrated easily with current operational implementations, thereby making this 
effort more likely to have a real impact. Finally, it should be noted that the reg- 
ularization constraint in field alignment is a weak constraint and the weights de- 
termine how strongly the constraints influence the flow field. The constraint in L 
is modeled as such because we expect the fluid flow to be smooth. From a reg- 
ularization point of view, there can be other choices ll27ll as well. The proposed 
method can be used for a variety of velocimetry applications including PIV, ve- 
locity from tracer-transport, and velocity from GOES and other satellite data, and 
an application of this is to advect rain-cells produced by a rainfall model, with 
realistic wind-forcing. 
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Figure 5: Deriving velocimetry information from satellite observations, Nexrad 
(top), GOES (bottom). See text for more information. 



